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On the Chebyshev approximation of a fnnction with two 

variables 

E. Scheiber* 


Abstract 

There is presented an approach to find an approximation polynomial of a func¬ 
tion with two variables based on the two dimensional discrete Fourier transform. 
The approximation polynomial is expressed through Chebyshev polynomials. An 
uniform convergence result is given. 
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1 Introduction 

The purpose of the paper is to present some aspects about the construction of 
an approximation polynomial for a function with two variables. The approximation 
polynomial is expressed through Chebyshev polynomials. 

Constructing an approximation polynomial of a function with the corresponding 
applications is the subject of the Chebfun software, presented in details in m, m- 
The Chebfun2 part of the software deals with the construction of an approximation 
polynomial of a function with two variables. According to 0, 0, to this end it is used 
a method based on Gaussian elimination as a low rank function approximation. 

In Chebfun the approximation polynomial of a function with one variable is obtained 
using an one dimensional discrete Fourier transform. The approach of this paper will 
use a two dimensional discrete Fourier transform. 

In spectral methods the Chebyshev polynomials are often used. The same form of 
the approximation polynomial are used in 0 , too. 

After recalling of some formulas on the Fourier series for a function with two vari¬ 
ables and on the two dimensional discrete Fourier transform there is presented an 
algorithm to obtain an approximation polynomial of a function with two variables and 
a convergence result. A Lagrange type interpolation problem for a function with two 
variables is studied. Two applications are mentioned; a numerical integration formula 
on a rectangle and a numerical computation of the partial derivatives. 
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2 Two dimensional Fourier series 


Let / : —)• i? be a continuous periodical function in each variable with the period 

27r. The Fourier series attached to the function is [8], t.3 

CX> 

f{x,y) ~ {an,m cos nx cos my + bn,m cos nx sin my+ 

ri,,m=0 


+Cn,m sin nx cos my 4 
with the coefficients given by 

aofi = ^ JJ^f{x,y)dxdy i 

anfi = ^ fffi f(x, y) cos nxdxdy i 

ao,m = ^ ffn f(x, y) cos mydxdy i 

ho,m = ^ JJq fix, y) sin mydxdy i 

Cn,o = ^ ffQf(x,y)sinnxdxdy 

where fl = [0, 27r]^. 

The complex form of the Fourier series i 

Z] 

n.mGZ 

with 

7 o,o = “0,0 


- dn,m sin nx sin my) 

^n,m = Unfix,!/) cos nx cos mydxdy 
'^n,m = ^ fix, y) cos nx sin mydxdy 
'in,m. = Unfix,!/) sin nx cos mydxdy 
in,m = U Unfix’d') sin nx sin mydxdy 

s 

^inx-\-imy 


UnU — 2(^^;0 
'T0,m — 2(^0,m 
Un,m — 4(^n,m 

Un^—m — “1“ ^n,m 


U — n^O — 2(^?^,0 “ 1 “ ^^ 72 , 0 ) 

To,—m — 2 {^ 0 ,m H“ ^^0,772) 

T—71,771 — 4(0^71,777 ^^71,777 “1“ i ^ n , 7 n ^77,777) 

U—n^—m — 4(0'?7,777 “1“ ^^77,771 “1“ i^n,m dn^m 


or 

7m,n = ^ y^/(x,y)e"*”"'"*™^dxdy, Mn,m^TL. (1) 

If the function is even in any variable then the hn,m,Cn,m,dn,m coefficients are all 
zero. 
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We shall suppose that the convergence conditions of the Fourier series to f{x,y) 
are fulfilled (the function has bounded first order partial derivatives in Q and in a 
neighborhood of {x,y) there exists or which is continuous in {x,y), cf. [8], 

t.3, 697). 

3 Two dimensional discrete Fourier transform 


Let be the inhnite matrix {xkj)k,j&z with the periodicity properties Xk+pj = Xkj, Xkj+q 
Xkj, yk,j G Z. The discrete Fourier transform construct another infinite matrix {yr,s)r,sez 
with an analog periodicity properties defined by 


p-l q-l 

yr,s = ^ ^ XkJ 
k=0j=0 


_• 27rfcr • 2Tr js 

e p e 1 . 


for r G {0,1,... — 1} and s G {0,1,... , g — 1}. 

The complexity to compute the pq numbers with the discrete fast Fourier transform 
algorithm is pq log 2 pq. 

As an application, if the Fourier series coefficients m are computed using the trape¬ 
zoidal rule for each of the iterated integrals then: 


ln,m = ^ I y)e—^ e—" fix, y)e-^^^dy"j dx 


^ J_ 

^4vr2yo 

Ve ^ / fix, -) 


27r^ 27rT -im^ 


ft 


1 dx = 


^ .=0 

j-m^dx p 


q 

p-l q-l 


1 ^TtT -in^ -im^ 

-,-)e p e <1 . 

k=0j=0 ^ ^ 


Thus the Fourier coefficients (7n,m) may be computed applying the discrete Fourier 

transform to f/(^, 

\ P ^ ) fcs{0,l,...,p—l},j£{0,l,...,g—1} 


4 The Chebyshev series 

Considering a continuous two real variables function /(x,y), x,y G [—1,1] there is 
attached the Chebyshev series 


CXD 

fix,y)^ ^ an,mTnix)Tmiy) (2) 

n,m=0 
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where 


ao,o - ^ fulfil 




dxdy 


a, 


_ 2 fi fi fi^,y)Tn(x) 


_ 2 

1,0 - J-iJ-i 


'JX-x^^X-y^ 


dxdy 


2 rr rr 

ao,m - J-iJ-i 


2 rl rl f{x,y)Tm{y) 


VX-x-‘^X-y'^ 

Changing x = cost,y = coss, the coefficient an,m will be 

4 


dxdy a.,. = ^^S^dxdy 


O^n.m. — 


n TT 

f (cos t, cos s) cos nt cos msdtds = 
~o /(cost, COS s) COS nt COS msdtds. 

J Jn 


( 3 ) 


Analogous formulas may be obtained for ao, 0 ) On.o and ao,™, too. Thus the coefficients 
of the Chebyshev series are the coefficients of the Fourier coefficients of the function 
(^(t, s) =/(cost, cos s). 

If the function / has second order derivatives then the Fourier series attached to (p 
converges to p and consequently 


f{x,y) = ^ an,mTn{x)Tm{y), x,yG[-l,l]. 

n,m=0 


(4) 


The polynomial 

n m 

fn,m{x,y) = '^'^ak,jTk{x)Tj{y) 
fc=0 j=0 

is called the Chebyshev approximation polynomial of the function /(x, y) in the square 

[-1,1]^. 

The parameters n,m are determined adaptively to satisfy the inequalities \akj\ < 
tol{= 10“^®, machine precision), for k > n and j > m. This is the goal of the algo¬ 
rithm [TJ The coefficients whose absolute value are less then tol are eliminated and the 
remained coefficients are stored as a sparse matrix. 

The fn,m{x,y) polynomial may be obtained with the least square method as the 
solution of the optimization problem 


mm / / 

J-xJ-X 


1 

^/l — — y'^ 



EE Xk,jTk{x)Tj{y) 

k=0 j=0 


2 

dxdy. 


Due to the Parseval equality 



n=0 m=0 n,m=l 



VT 


f{x,y) 

- x2^1 - y2 


dxdy 
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the quality of the approximation polynomial may be evaluated by 


J_ fix,y) 

7r2 vi-xVi-y^ 



Algorithm 1 Algorithm to compute the Chebyshev approximation polynomial 
1: procedure CHEBFUN2(f) 

2: n •(— 8 

3: tol ^ 10-15 

4: sw ■(— trrte 

5: while sw do > The approximation polynomial is determined adaptively 

6: m 2n 

7: x,y ^ cos —, fc = 0 : m — 1 

8: z^f{x,y) 

9: g FFT (z) /m? 

10: a ^ 45Rg'(l : n, 1 : n) 

11: a(l,l)^a(l,l))/4 

12: a(l, 2 : n) •(— a(l, 2 : n)j2 

13: a(2 : n, 1) •(— a(2 : n, l)/2 

14: if |a(i — 1 : i,l : n)\ < tol &Z, |o(l : n,i — 1 : i)\ < tol then 

15: sw •(— false 

16: else 

17: n ^ 2n 

18: end if 

19: end while 

20: for i = 1 : n do i> Removal of negligible coefficients 

21: for j = 1 : n do 

22: if \aii,j)\ < tol then 

23: ^ 0 

24: end if 

25: end for 

26: end for 

27: return a 

28: end procedure 


The value of the polynomial fn,m in a point (x, y) may be computed adapting the 
Clenshaw algorithm, [9], but we find that the evaluation of the expression fn,m{x,y) = 
Vf{x)An,mVmiy), where 14(s) = (ro('S),ri(s), • • • and A = {ak,j)k= 0 :n,j=o-.m is 

more efficient within a matrix oriented software. V' denotes the transpose of the vector 
V. The complexity order of both algorithms is 0{nm). 
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5 The Chebyshev series of partial derivatives 


We assume that the function f{x,y) has first order continuous partial derivatives 
and the series (jl]), there is required to find the coefficients {bn,m)n,meN such that 


Yl KmTn{x)TM- 


n,m=0 


Using the equalities 
and 


T[{x) = To{x), r'(x)=ri(x) 
1 [T'^+i{x) Tl^-i{x)' 


n + 1 


may be written as 


dfix,y) 

dx 


n — 1 


OO /CO 


= T„(x), n > 1. 


'y ^ ( y ^ bn^rnTnjx) j Tm{y) 


m=0 \n=0 


A ^ hrnnix) ^bkm TLlix) U_i{x) 

Y 1 KmTiix) + +Y 

k=2 


m=0 

OO 


k-1 


) Tm{y) = 


Y (^ 0 -^ 


m=0 


^)T{(x) + f; ^{bk-i,m - bk+i,Mx)] TM = 

k=2 ^ / 


CO /CO 


dfix,y) ^ , 

- 5 - = > Cn,mTn(x)Tm[y) 

dy ^ 


n,m=0 


satisfies the relations 


^ri,2 

Cn,0 9 

2j 


— Q^n,l 


(^nj —1 C77,j-|_i) — 


a 




J>2, 


m E N. 


( 6 ) 


— Y^ ( I '^rn{y)- 

m=0 \A:=1 / 

Identifying the coefficients of Tj^ix) there is obtained the linear algebraic system 

r ^2,m _ ^ 

^0,m 2 — ^l,m 

l,m ^ ^ 2, ?Tl G N. 

Summing the above equalities for A: = n + 1, n + 3, n + 5,... it results 

bn,7n — “t“ “t“ 3)(ln-\-3,m H“ “t“ 5)tt7T,-(_5^j7T, + . . . V U, G N. 

In the same way it is deduced that the coefficients of the series 


(7) 


( 8 ) 
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Let be 


(9) 


_ dfix,y) . 

dxdy ~ ^ CLn,m-Ln[.X}lm{y}- 


dydx 


n,m=0 


Because applying ([8]) it results 


1 dx 


Q ■ (^fcj —1 *^fc,j-|-l) ^k,j 

and consequently, for k,j > 1, 

i^k—i,j—i ^fc—ij'+i 1 “1“ ^fc+ij+i) (^9) 

Denoting = (ifc_ij_i—dfc+ij-i+dfc+ij+i), from the Parseval equality 
corresponding to ([9]) it results that 


EEAi<i6M?,i. 

k=2 j=2 


where Mi,i > max{\^^{x,y)\ : x,y e [-1,1]}. 

6 The convergence of the Chebyshev series 


Using the techniques presented in [7] and [1], in some hypotheses it may be proven 
that the convergence in Q is uniform in [—1,1]^, for n,m^ oo. First we state 

Theorem 6.1 If the function f has second order continuous derivatives then 


2M2 0 8M2 0 

1^71,o| — 7 TTT" and E 7 ryT"’ n ^ 1, 


(n — 1)^ 


7 r(n — 1 )^' 


i„, I ^ 2Mo,2 1 ^ I ^ 8Mo,2 _ ^ T 

|Q^0,m| — / T\o CLTid 0^1 m S / i\9’ TTl > 1, 

[m — ly 7r(m — ly 

where M 2 ,o > max{|0(x,y)| : x,y S [-1,1]}, Mo ,2 > max{|0(x,?/)| 

[- 1 , 1 ]}- 


( 11 ) 

( 12 ) 

x,y E 


Proof. The coefficient 0^,0 may be written as 

«nO = —7 / ( [ /(cos t, cos s) cos nsdt ) ds. 

^ /o V/o / 
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Two partial integrations are performed in the internal integral 


/•TT 

/ /(cost, cos s) cosntdt = 

Jo 


1 5^/, , /sin (n—lit sin(n+l)t\ , 

' - (cost, cos s) sin t 1 -^- V .. 1 dt. 


2n Jo 


n — 1 


It results that 


PIT ^ / 1 

/ /(cos t, cos s) cos ntdtl < —M 2 o ( - 

Jo 2n ’ V’^- 


n + 1 


_l_ ^ ^ < ^^ 2,0 


— 1 n + ly (n — 1) 


and consequently |an,o| < • 

Using (fT3]l in 0^,1 = fj (JJ^/(cos t, cos s) cos ntdt) cossds it results 




4 M2,0 


The proof of (I12p is similar. 


PTT 

/ I cos s|ds = 

Jo 


8 M 


2,0 


7r(n — 1)2 


(13) 


Theorem 6.2 If the function f has second order continuous partial derivatives then 
lim„, /n,m = / uniformly in [-1,1]^. 


Proof. From 

00 n m 

f{x,y)= ^ an,mTn{x)Tm{y) and fn,m{x,y) = '^'^ak,jTk{x)Tj{y) 
n,m =0 k =0 j =0 

it results 

n 00 00 00 

fix,y) - fn,m{x,y) = EE ak,jTk{x)Tj{y) + E E ak,jTk{x)Tj{y). 

k= 0 j=m+l k=n+lj =0 


Then 


n 00 00 00 


\f{x,y) - fn,m 

VI 

E: 

E 1 

^ ^ “fcj 



II 

0 

II 

=m+l 

/c^n+l 

00 

= E 

+ 

o' 

00 

E 


n cxD 

j=m+l 

j 



^=2 j=m+l 

00 

+ E 

Q^fc,o + 

00 

E 

afc,i 

00 00 

1 + Z] 


k=n+l fc=n+l fc=n+l j =2 
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and using the Cauchy-Buniakowsky-Schwarz inequality it follows 


{f{x,y) - fn,mi^,y)f < 


n oo 


< 


6 X] + X] + X] X] + 


\j=m+l 


L j=m+l 


k=2 j=m+l 


oo CO 




\k=n-\-l 


\k=n-\-l 


fc=n+l j=2 


The following inequality holds ^ ^ = ^• 

Using the results of Theorem 16.11 the first, second, fourth and the fifth expression 
are increased by 


CO 

I“Ob I - 

j=m+l 

oo 

Y 1 I “hi I < 

j=m+l 

CO 

i«fe,oi < 

k=n+l 

CO 

/c=n+l 


2Mo,2 Y.~2 < 

j=m 

8Mo,2 1 

TT p 

j=m 

OO ^ 

2M2,o Z] ^ < 

k=n 

8M2fi 1 

TT k'^ 

k=n 


2Mo^2 
m — 1 


8-^0,2 4Mo,2 

7r(m — 1) m — 1 

2M2,o 
n — 1 


8^2,0 ^ 4M2,o 
7r(n — 1) n — 1 


For the third and sixth expression we shell use (IlOp and then the Cauchy-Buniakowsky- 
Schwarz’ inequality 


n CO 


n oo 


E E K, i = E E 


.k=2 j=m+l 


k=2 j=m+l 


I Afcj 

4kj 


^ / n CO \ ( ^ i\ 

- 16 ( ^ ^ ) ( ^ ^ ^ ) 

yfc=2j=m+l j \^k=2j=m+l j 


< 


”1 “ 1 


k=2 j=m-\-l 


6m 


and respectively 


CO CO 


OO CO 


E Ei“«i < E E 


./c=nH-l j=2 


.fc=n+l j=2 


|Afcj 

4fcj 


£vi. E 


\k=n-\-l 






1 \ vr^M? 


2 I 6n 


1,1 
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Consequently 


\f{x,y) - fn,m{x,y)\ < 
when m,n ^ oo. m 


20M2 20M2 vr^M^ 


0 , 


7 The Lagrange interpolation polynomial 

For any grids —1 < xq < xi < ... < < 1, —1 < yo < ?/i < • • • < l/m < 1 and any 

/ : [—1,1]^ —^ M the expression of the Lagrange interpolation polynomial is 

n m 

Ln,m{x,y) ='^'^f{xi,yj)lxi{x)lyj{y) 
i=0 j=0 


where 


Ixi(x) = 


X - Xk 


, and lyjix) = 


y-yi 


k=0,k^i ‘ 


Xi X]^ 

This polynomial satisfies the interpolation restrictions 

Ln,mixk,yi) = f{xk, yi), VA; G {0,1,..., n}, and V/ G {0,1,, m}. 

In the set of (n, m) degree polynomials there exists a unique interpolation polynomial. 

If Xi = cos i G {0,1,..., n} and yj = cos j G {0,1,... , m} then using the 
discrete orthogonality relations, [9], 


0 if p / y 

^^lkTp{xk)Tq{xk) = I f if p = g G {1,2,... ,n - 1} = nap6p^ 
k=o ( n if p = gG{0,n} 




where 


1 u 


, 2 if k£{0,n} _( I if iG{l,2,...,n-l} 

if A;G{l,2,...,n-l} ana | ^ iG{0,n} 

the Lagrange interpolation polynomial may be written as 

n m 

Ln,m{x,y) = CijTi{x)Tj [y), 

i=0 j=0 

where Qj = :;^ln,ilra,jTJk=o'llT=oln,klra,lf{xk,yi)Ti{xk)Tj{yi). 

This polynomial will be called the Lagrange-Chebyshev interpolating polynomial. 
As in [4], the following statements occur: 
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Theorem 7.1 (Aliasing of Chebyshev polynomials, W) For any n > 1 and 0 < m < n 
the polynomials Tm,Tn±m,T 2 n±m, ■ ■ ■ take the same values on the grid (cos ^)o<fc<n- 

Theorem 7.2 (Aliasing formula of Chebyshev eoeffieients) Let 

OO OD 

f{x,y) = J2'^^hjTi{x)Tj{y) 

i=0 j=0 

and let Ln,m{x,y) = Yll(=o^'f=o^i,j'^iix)Tj{y) he its Lagrange-Chebyshev interpolant. 
Then 

CO CO 

^ ^ ay.2pn+i,2qm+j “1“ ^ ^ (^2pn—i,2qm—j (1^) 

p,q=0 P)9=l 


Proof. Supposing that {cij)o<i<m,o<j<m are given by ([l5]) and (p{x,y) = 
Er=oEiloCijTi(x)rj(y). For any {k,l) G {0,1,..., n} x {0,1,..., m} 

n m 

T{xk,yi) = EE=m Ti{xk)Tj{yi) = 

2=0 j = 0 


n m / OO 


= EE E C^2pn-\-i,2qm+j ^ ^ Oi2pn—i,2qm—j 1 'Li{,Xk')Tj(^yif 

1=0 j=0 \n,q=0 Pi9=l / 

It is observed that when the indexes i,j,p, q go through their values then (2pn±i, 2qm± 
j) go through N X N and any two pairs are distinct. Thus, with Theorem 17.11 

OO OO 

T{xk,yi) = EE“.. tTs{xk)Tt{yi) = f{xk,yi). 

s=0 t=0 

The unicity of the interpolating polynomial in the set of (n, m) degree polynomials 
implies Ln^m = T- ■ 

A consequence of (fTSl) is a relation between Ln^mix, y) and the approximation poly¬ 
nomial fn,ra{x,y) = Er=0 E^Lo (?/) ^ 

n m 

Ln,m{x,y) = ^ ^ ^ ^^ c-i,jTdi{x)Tj(y') 
i=0 j=0 


EE + E Oii,2qm±j T ^ ^ C^2pn±i,j T ^ ^ 0!2pn±i,2qm±j 1 'ddi{x)Tj(gj) — 


i=0 j=0 


9=1 


p=l 


p,q=l 

OO m 


/»,™{^.rf + E E Oii,jTi{x)Tp.{y) + E E aijTy^{x)Tj{y)+ 

2=0 j=m-\-l 2 = 71+1 j = 0 
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Then 


CO CO 

+ 0'ijTiy.{x)Tf^^{y). 

2=n+l j=m+l 


y)\^Yl Y1 E E \aij\. ( 16 ) 

i=0 j=m+l i=n-\-l ^=0 

Know, it can prove the uniform convergence of the Lagrange-Chebyshev interpola¬ 
tion polynomials: 

Theorem 7.3 If the function f has second order continuous partial derivatives then 
\\mn,m^oo Ln,m = f Uniformly in [- 1 , 1 ]^. 

Proof. Using (fTT)l and (fT6l) we obtain 

Ifi^iU) -bn,m(3:,y)| E \f{x,y) fn,m\ T \fn,m(,X,y') Ln,m{,X,y)\ E 


n CO 


CO CO 


£2 E E !“«!+ E E 

\i=0j=m+l i=n-\-l j=0 

The rest of the proof follows the proof of Theorem [ 


8 Applications 

1. Integrating fn,mix,y) on O there is obtained 


f{x,y)dxdy : 


fn,m{x,y)dxdy = A ^ ^ 


CXr 


k=0, evenj=0, even 




;5T' (17) 


2. Computation of the first order partial derivatives. Practically, knowing the 
Chebyshev approximation polynomial fn,mix, V) = Y2=oYlJ=o(^k,jTk{x)Tj{y), and 
with the assumption that ak,j ~ 0 for A: > n > 4, and for any j E {0,1,...,m} 
the hrst n equations of the system ([7]) will be 


with the solution 

bn—l,j 

bn—2.j 

bk,j 

bo,j 


Ki - ¥ 

= “1,1 


^(bk-ij — bk+i,j) 

— 

/c E {2,3,..., n — 2} 

1 

2(n-l) 

— O^n—lJ 


1 K 

2n "n-l.l 

= Oin,j 


= 2,TlQinJ 
— 2(?1 
— 2(fc -|- 

= “1,1 + 

+ bk+2,j, 

A: E {n — 3, n — 4,..., 2,1} 
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Then ^ 

The partial derivative may be computed similarly. 

Due to the truncation of the Chebyshev series the numerical result is influenced by 
the truncation error as well as by rounding errors. The automatic differentiation [3] 
is a method which eliminates the truncation error but it requires a specific computa¬ 
tional environment related to the definition of the elementary functions (e.g. apache 
commons-mathS v. 3.4). 

9 Examples 

Using a Scilab implementation the following results are obtained 
1. f{x,y) = cos XT/ [2], Ch. 11, p. 2. The matrix of the coefficients is 


0.880725579 

0. 

- 0.117388011 

0. 

0.001873213 

0. 

0. 

0. 

0. 

0. 

-0.117388011 

0. 

- 0.114883808 

0. 

0.002484444 

0. 

0. 

0. 

0. 

0. 

0.001873213 

0. 

0.002484444 

0. 

0.000603385 


The value of the indicator given by ([5]) is 3.97247 • 10 

On an equidistant grid of size 50 x 50 in [0,1]^ the maximum absolute error is 
0.000082141. 

The integral given by (fTTll is 3.784330902, while Mathematica gives 4SinIntegral[l] : 
3.78433228147. 

2. g{x) = cos lOxy^ + e~^^ [2], Ch. 12, p. 6. The size of matrix of coefficients is 
33 X 43. 

The value of the indicator given by ([5]) is 0. 

On an equidistant grid of size 50 x 50 in [0,1]^ the maximum absolute error is 
2.98594 • 10-i3_ 

The integral given by (11711 is 4.590369905, which is equal to that given by Math¬ 
ematica. 

10 Conclusions 

There is presented an alternative to the Gaussian elimination method used in Cheh- 
fun software in order to construct an approximation polynomial of a function with two 
variables. 

Because the discrete Fourier transform is a common tool for the usual mathematical 
softwares, this approach has a relative simple implementation, but as a drawback, if 
the tolerance is the machine precision then it may require a large amount of memory. 
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